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We consider pair production and decay of fundamental unstable particles in the framework of a 
modified perturbation theory (MPT), which treats resonant contributions of unstable particles in 
the sense of distributions. The cross-section of the process is calculated within the NNLO of the 
MPT in a model that admits exact solution. Universal massless-particles contributions are taken 
into consideration. The calculations are carried out by means of FORTRAN code with double 
precision, which ensures per mille accuracy of computations. A comparison of the outcomes with 
the exact solution reveals an excellent convergence of the MPT series at the energies close to and 
above the maximum of the cross-section. Near the maximum of the cross-section the discrepancy 
of the NNLO approximation makes up a few per mille. 
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A description of processes of production and decays of fundamental unstable particles for 
colliders subsequent to LHC should be made generally with the NNLO accuracy. This implies 
that the calculations must provide gauge cancellations and unitarity, and enough high accuracy 
of computation of resonant contributions. The existing methods, such as the DPA successfully 
applied at LEP2 ^ or the complex-mass scheme (CMS) [Q] intended mainly to ILC [gj, provide 
only the NLO precision of the description. The pinch-technique method [Q], another for a long 
time developed approach, in principle can provide the NNLO precision, but for the maintenance of 
the gauge cancellations it requires a lot of calculations of extra contributions that become apparent 
formally at the next level of the precision, which is impractical [Q]. So alternative approaches are 
required. A modified perturbation theory (MPT) [Q, [7|, ^ is such an approach. Its main feature 
is the direct expansion of the probability instead of amplitude in powers of the coupling constant 
with the aid of distribution-theory methods. As the expansion is made in powers of the coupling 
constant and the object to be expanded is gauge invariant, the gauge cancellations in the MPT 
should be automatically maintained. However, the accuracy of the description remains unknown. 
To clear up this question numerical simulations are required. In this report I present results of the 
appropriate calculations and, simultaneously, do a brief introduction to the MPT method. 

Actually I will discuss the case of pair production of unstable particles. The most crucial in 
this case are the double-resonant contributions. The corresponding total cross-section, for example 
in e + e~ annihilation, has the form of a convolution of the hard-scattering cross-section with the 
flux function, 



(In fact the angular distributions may be described in the MPT, as well, but I do not discuss this 
option below.) The hard-scattering cross-section has the form of an integral over the virtualities of 
unstable particles, 



a(s)= / / dsids 2 d{^fs-y/J[-^)^X{s,s i ,s 2 ) <Z>(s;s u s 2 ) p(si)p(s 2 ) ■ (2) 



Here the first two multipliers in the integrand are kinematic factors, p(s,) are Breit-Wigner (BW) 
factors, (\ + 8 c ) stands for soft massless-particles contributions, and function <J> is the rest of the 
amplitude squared. Generally <I> corresponds to one-particle irreducible contributions, and hence <I> 
has no singularities on the mass-shell of unstable particles. On the contrary, the kinematic factors 
have singularities due to the theta function and the square root of the kinematic function A. The 
BW factors if to naively expand them in powers of the coupling constant a generate nonintegrable 
singularities. Actually this constitutes a serious problem since integrals in @ become senseless. 

The singularities, nevertheless, become integrable if to expand the BW factors in the sense 
of distributions. In this case the expansion of a separately taken BW factor is beginning with the 
8-f unction which corresponds to the narrow- width approximation. The contributions of the naive 
Taylor expansion are supplied with the principal-value prescription for the poles. The nontrivial 
contributions are the delta-function and its derivatives with coefficients c n , which are polynomials 
in a determined by the self-energy of the unstable particle [^]. Within the NNLO, the expansion 
has the form 
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Here M is the renormalized mass, To is the Born width, L(s) is the self-energy of the unstable parti- 
cle. Coefficients c n within the NNLO include 3-loop self-energy contributions and their derivatives, 
determined on-shell. The structure of the contributions is such that in the OMS-type schemes of 
the UV renormalization the real self-energy contributions enter into the coefficients either without 
the derivatives or with the first derivative only. This implies that the relevant real self-energy con- 
tributions are determined by the renormalization conditions. Unfortunately, the conventional OMS 
scheme is not convenient, since it does not maintain the gauge independence of the renormalized 
masses of unstable particles. So it is reasonable to proceed to the particular version of the OMS 



scheme, namely to the OMS or pole scheme []9[ |10p, where the renormalized masses of unstable 



particles by definition coincide with the real parts of the poles of the propagators and thereby co- 
incide with the observable masses. The coefficients c n in this scheme are determined as follows 



7 2 , 2 

fo = - a -r + a 
h 
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Cl =0, c 2 = -a 2 l\. (4) 



Here I k = Iml^M 2 ), L = a'L l + a 2 L 2 + a 3 I 3 , and l[ = ImIj(M 2 ). Simultaneously in the OMS 
scheme the Im£(M 2 ) coincides with the imaginary part of the pole of the propagator. This allows 
one to connect 4 order-by-order with the width of the unstable particle via the unitarity relations 
ah =MT , a 2 I 2 =MaT u and a?h = Ma 2 T 2 + T 3 /(SM) [§. 

In fact, expansion (||) has sense only if the weight in the integral is a regular enough function. 
In our case, however, the kinematic factor is not regular, and this leads to the divergence of integrals 
(||) after the substitution of the expansions. At first glance this closes a possibility of application 
of the expansions for the BW factors. Nevertheless, the kinematic factor may be analytically reg- 
ularized via the substitution [X(s,si ,S2)] 1 ^ 2 —> [X(s,si,S2)] v . With large enough A this imparts 
enough smoothness to the weight, and the singular integrals become integrable. Moreover, after 
the calculation of the integrals and removing the regularization the outcomes remain finite and the 
expansion remains asymptotic [B]. In principle, this salvages the applicability of the approach. 

However, the singular integrals must be analytically calculated. The scheme of their calcu- 
lation is as follows. At first one should proceed to dimensionless energy variables x, Xj (i = 1,2) 
counted off from thresholds, ^fs = 2M(1 +x/4), yfsi = M(l +x,/2). Formula (§) for the hard- 
scattering cross-section (with the analytically regularized kinematic factor) then takes the form 

o(x) = dxidx 2 (x-x l -x 2 )\^>(x;xi,x 2 ) p(xi)p(x 2 ) . (5) 



Here (x— x\ — x 2 ) v + = 6{x—x\ —x 2 ){x—x\ —x 2 ) v and tilde marks the dimensionless functions. (For 
convenience, factor (1+5 C ) is included into the definition of the test function <!>.) Further, I sub- 
stitute asymptotic expansions for p(x,) and consider at every rij in the product of the expansions 
the contributions of the PVxJ ni and 8 ( - n ^ 1 \xi). Simultaneously, in each case, I represent the test 
function in the form of double Taylor expansion truncated at x["' l ' and x^ 2 l ' plus a remainder, 

^ (i|-ln 2 -l Ji x k2 ____ 

$(jc;jci,jc 2 ) = L E 7^ 7^7 <S> (k, ' k2) (x;0,0) + AO(x;x 1 ,x 2 ) . (6) 

*i =0*2=0*! ! k2 ' 
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The higher powers of xi and x% in the Taylor expansion will zero the 8^ n ~^(xi) and cancel the 
PVxJ m . The remainder A<I> is determined as the difference between <1> and the Taylor expansion. 
In fact A<i> is to be further expanded with respect to separately x\ and x%, but for brevity I do 
not consider this procedure explicitly (see details in [Jsj] ) - I mention only that the final remainder 
produces a regular contribution to the integrand in formula (||), and the integrals of this contribution 
can be numerically calculated. At the same time, the contributions of the Taylor expansions are 
singular, and the integrals of them are analytically calculable. After the analytic integrating and 
putting v — ^ 1/2, a sum of regular and singular contributions appears with the singular contributions 
being products of regular factors and the power distributions of the type x 5 / 2 " with integer n. 

The convolution integral ([]]) of above contributions can be numerically calculated. In particu- 
lar, the integral of singular power distributions may be calculated by means of the formula 




where (j) is a weight and A is a positive integer such that —N— 1 < Re v < —N. The test function 
<1> at this stage is considered determined in the conventional perturbation theory. 

For carrying out the above mentioned calculations rather general FORTRAN code with the 
double precision is written. The calculation of regular integrals in this code is fulfilled by Simpson 
method. Numerous indeterminate forms of the type 0/0 that emerge in the integrand due to the 
difference structures are resolved through the introduction of linear patches. The patches diminish 
the errors that arise because of the loss of decimals near the indeterminacy points, and manifest 
themselves, in particular, in numerical instability. Of course, the patches generate the errors, too. 
However the latter errors are under control and are estimated as £ 2 <Pq/(3<^), where e is the size 
of the patch, (po and <pQ is the integrand and its second derivative at the center of the patch. This 
estimate has sense of a relative error of the calculation of integral inside the patch. The sizes of 
the patches are so determined that the singular components of indeterminate forms take a constant 
value on the boundaries of the patches. So, for instance, in the case of an indeterminate form of the 
type x" fx" the size e should be so chosen that e" would be a constant, say, 10~ N with some fixed A". 

Unfortunately, the errors caused by the loss of decimals near the indeterminacy points but 
outside of the patches, cannot be explicitly estimated. What is possible to formulate is only the 
dependence of these errors on the basic parameters. So, in the case of indeterminacy [f(x) — 
f(0)]/x the appropriate error has a behavior ~ lO~^ D ~ N ^fo/fQ, where D is the number of digits in 
the representation of real numbers (in our computations D = 15), and A" is the above mentioned 
parameter of smallness of e. In the case of [f(x) — /(0) — xf'(x)]/x 2 , the analogous estimate is 
~ 10~( D ~ N ' fo/ f$ , where N is the order of smallness of e 2 etc. 

Nevertheless, the total error of the computations can be estimated. The crucial point is that the 
errors because of the patches are increasing with increasing the sizes of the patches, while the errors 
because of the loss of decimals are decreasing. This implies that there should be an optimum size of 
the patches when the sum of all the errors is minimized. The point of the minimization must possess 
extremum properties, so that the result of the computation at this point must be stable with respect to 
varying the sizes of the patches. Moreover, at the extremum point, the sums of the errors of different 
kinds must be approximately equal each other (up to a coefficient of order one). So the order of 
the total error of the computation may be estimated by the order of the sum of errors because of 
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Figure 1: The total cross-section: the exact result in the model and the results of computations in MPT. 
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the patches. However, the latter sum can be numerically estimated. By this means the total error 
of the computations can be estimated. In particular, in the case of the model discussed below, the 
optimum size of the patches is found established by N « 8 (at D = 15). The corresponding estimate 
of the relative error of the computation of the NNLO approximation constitutes 10 . 

Now I proceed to the physical model in the framework of which I do computation. First I recall 
that my aim is to verify whether the MPT calculations are practicable in principle, and then my aim 



is to test the convergence properties of the MPT expansion. The existing experience [ ]11[ ] shows 
that the convergence properties are determined mainly by the BW factors, but not by particular 
form of the test function. So without loss of generality I can determine the test function <J> in the 
framework of a model. As such a model, I consider the Born approximation for the process e + e~ — > 
y,Z —>tt — >■ W + bW~b. Simultaneously, I consider the self-energies in the denominators of the top 
quark propagators with the 3-loop contributions [jjl]]. Since soft massless particles generally can 
make considerable contributions, I take into consideration also the universal soft massless-particles 
contributions. Actually they are collected in the flux function and in the Coulomb factor. The flux 
function, I consider in the leading-log approximation. In the Coulomb factor determined mainly by 
gluons, the multi-gluon contributions, generally, are important. However, at least at a distance from 
the threshold a qualitative picture may well be simulated by the one-gluon contributions. So, in the 
framework of a model I consider Coulomb factor in the one-gluon approximation. In addition, I 



imply the standard resummation in the Coulomb factor Q12|], which does not affect the BW factors. 

The outcomes of the computations are presented in the Figure and in the Table. In the Figure in 
the panel (a) the thick curve shows the behavior of the total cross-section in the model. The dotted, 
dashed, and continuous thin curves show the results of the MPT computations in the LO, NLO, 
and NNLO approximations, respectively. It is worth noting that the NNLO result almost coincides 
with the exact result in the model. The distinction is visible in the panel (b) where the percentages 
with respect to the exact result are presented. In the Table the outcomes are represented in the 
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numerical form at the characteristic energies planned at the ILC [Qj. In the last column the numbers 
in parenthesis represent the uncertainties in the last digit. In other columns the uncertainties are 
omitted as they appear beyond the precision of the presentation of data. 

In conclusion, first of all the above results show in practice the existence of the MPT expansion 
in the case of pair production and decay of fundamental unstable particles. Secondly, the NLO and 
NNLO approximations in the MPT have very stable behavior at the energies near the maximum of 
the cross-section and at higher energies (and this result to a large extent is model-independent [11]). 
Thirdly, at the ILC energies the NNLO approximation in the MPT in the case of the top-quark pair 
production provides approximately a half-percent accuracy of the description of the cross-section. 
In fact this is what is needed at the ILC. However, the higher precision of the MPT description, 
when needed, may be achieved by the proceeding to the NNNLO, or by the using of the NNLO for 
the calculation of the loop contributions only — on the analogy of actual practice of application of 
DPA [@]. The higher precision of the actual computations in the NNLO of MPT may be achieved 
at the proceeding to the representation of real numbers with higher precision (higher than with the 
double precision). 

The author is grateful to A.L.Kataev for encouragement and suggestion to make a report at the 
ACAT Workshop. 
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